Efficient Bayesian Hierarchical Modeling of Spatial Precipitation Extremes

An efficient Bayesian hierarchical model for spatial extremes on a large domain is proposed. In the data layer a Gaussian elliptical copula having generalized extreme value (GEV) marginals is applied. Spatial dependence in the GEV parameters are captured with a latent spatial regression. Using a composite likelihood approach and a method for incorporating stations with missing data, we are able to efficiently incorporate a large precipitation dataset. The model is demonstrated by application to seasonal precipitation extremes at approximately 2800 stations covering the western United States, -125E – -100E longitude and 30N to 50N latitude. The hierarchical model provides parameters on a 1/8th degree grid and consequently maps of return levels and associated uncertainty for each season. The model results indicate that return levels vary coherently both spatially and across seasons, providing valuable information about the space-time variations of risk of extreme precipitation in the western US, helpful for infrastructure planning. Engineering design of infrastructure such as flood protection, dams, etc. and management of water supply and flood control require robust estimates of return levels and associated errors of precipitation extremes. Spatial modeling of precipitation extremes not only can capture spatial dependence between stations but also reduce the overall uncertainty in at-site return level estimates by borrowing strength across spatial locations [1]. Hierarchical Bayesian modeling of extremes precipitation was first introduced by [1] and since has been widely discussed in the literature [2], [3]. Hierarchical modeling is an alternative to regional frequency analysis providing gridded or pointwise estimates of return levels within a study region [4]. Bayesian hierarchical models for spatial extremes have typically been limited to small geographic regions that include on the order 100 stations covering areas on the order of 100,000 km2. Large geographic regions with many stations present a computational challenge for hierarchical Bayesian models, especially when computing the likelihood of Gaussian processes (GPs), which for Corresponding author: C Bracken, University of Colorado, Boulder CO [email protected] 1 Bureau of Reclamation 2 Cooperative Institute for Research in Environmental Sciences 3 n data points, requires inverting an n × n matrix, an O(n3) operation. Several approaches exist for speeding up GP likelihood computations such as low-rank approximations [5] in which the GP is approximated at a small number of knots and composite likelihood methods [6] where the likelihood computation is broken into groups containing a small number of stations. The use of a composite likelihood approach is explored here because we not only wish to estimate covariance parameters but to also produce maps of return levels with small credible intervals. Some attempts have been made to model extremes in large regions and large datasets in a Bayesian hierarchical context. [7] use a hierarchical max-stable model with climate model output in the east coast to examine spatially varying GEV parameters, with a max stable process for the data dependence level. [8] model gridded precipitation data over the entire US, for annual maxima at a 5x5 degree resolution (43 grid cells) and copula for data dependence, incorporating spatial dependence directly in a spatial model on the data, not parameters. [9] and [10] model over 1000 grid cells of climate model output using spatial autoregressive models which take advantage of data on a regular lattice to simplify computations. The research contributions of this study are as follows. A Bayesian hierarchical model is proposed which is capable of incorporating thousands of observation locations by utilizing a composite likelihood method. The GEV shape parameter is modeled spatially in order to capture the detailed behavior of extremes in the western US. In addition the model is capable of incorporating stations with missing data with little additional computational overhead. The model is applied to observed precipitation extremes in each season, providing estimated seasonal return levels for the western US.